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A  NUMERICAL  SCHEME  FOR  PREDICTING  THE  LOCATION 
OF  TIDALLY-GENERATEO  FRONTS  IN  SHALLOW  WATER 


Alan  J.  Elliott 


ABSTRACT  , 

~^y  .  /X  +V  yyi* 

A  hydrodynamic  manerica l Jmode l  is^used  to  compute  tidal  currents  and 
calculate  the  parameter  (log^\({0/h)  ,  where  u  is  the  maximum  local 

tidal  speed  and  h  is  tne~depth.  A  recent  theory  has  shown  that 
fronts ,  marking  the  boundary  between  stratified  and  isothermal  water, 
can  be  expected  to  form  where  this  parameter  takes  a  critical  value  of 
about  2.2.  Thus,  if  the  amplitude  and  phase  of  the  tidal  elevations  are 
known  around  a  region  of  interest  then  the  likely  location  of  fronts  can 
be  predicted  if  the  bottom  topography  is  known.  As  an  example  of  the 
method  the  scheme  is  applied  to  the  Southwestern  Approaches  to  the 
English  Channel ,  and  good  agreement  is  obtained  for  the  predicted 
location  of  a  thermal  front  that  has  been  observed  at  the  western 
entrance  during  the  summer  months.  In  particular ,  the  model  shows  that 
even  eastward  winds  of  Beaufort  force  8,  which  occur  for  at  least  10% 
of  the  time  in  this  region,  are  unlikely  to  influence  the  location  of 
the  front.  The  computer  program  and  instructions  for  its  use  are  given 
in  an  appendix. 


INTRODUCTION 

During  recent  years  there  has  been  a  growing  body  of  evidence  to  suggest 
that  thermal  fronts  can  be  generated  In  shallow  water  by  the  bottom¬ 
generated  turbulence  that  Is  associated  with  tidal  currents. 

The  first  observations  of  this  effect  were  made  by  Simpson  [ij  who 
observed  a  stratified  region,  persistent  throughout  the  summer  months. 

In  the  northwestern  part  of  the  Irish  Sea.  This  region  coincided  with  a 
zone  characterized  for  Its  weak  tidal  currents,  the  remainder  of  the 
Irish  Sea  being  noted  for  Its  strong  tidal  currents  and  vertically  well- 
mixed  water.  In  a  further  Investigation  [2]  it  was  shown  that  the  rate 
at  which  work  Is  done  by  the  tidal  currents  against  the  bottom  friction 
Is  proportional  to  us,  where  u  Is  the  tidal  velocity,  and  that  the 
criteria  on  whether  the  bottom-generated  turbulence  will  be  sufficiently 
strong  to  completely  mix  the  water  column  vertically  against  buoyancy 
forces  will  depend  on  the  ratio  of  u*/h,  where  h  Is  the  local  bottom 
depth.  By  using  a  numerical  model  to  predict  the  tidal  currents  through¬ 
out  the  Irish  Sea,  and  by  comparing  the  contoured  values  of  logiQ(us/h) 

with  the  locations  of  known  fronts,  Simpson  and  Hunter  showed  that  the 
critical  value  for  the  frontal  location  was  around  2.2  (In  cgs  units). 


1 


This  work  was  extended  by  Feewnhead-  /i/  who  used  mean  tidal  charts  to 
construct  the  contours  of  the  stratification  parameter  In  the  waters 
surrounding  the  coast  of  the  British  Isles.  Similar  calculations  of 
this  kind  have  also  been  made  by  using  tidal  charts  for  conditions  of 
mean  spring  tides  and  comparing  the  results  with  the  frontal  locations 
determined  from  satellite  Images  [4.  57.  In  a  detailed  numerical  study 
of  a  front  observed  In  the  southern  Irish  Sea,  James  [6]  showed  the 
seasonal  development  of  the  frontal  system  and  demonstrated  that  a  front 
that  forms  during  the  neap  part  of  the  tidal  cycle  may  not  be  destroyed 
during  the  following  spring  tide.  More  recently,  Plngree  and  Griffiths 
17]  have  used  a  high  resolution  numerical  model  to  calculate  the  strati¬ 
fication  parameter  on  the  continental  shelf  around  the  British  Isles, 
and  good  agreement  was  obtained  between  the  predicted  and  observed 
frontal  locations. 

The  purpose  of  this  memorandum  Is  to  present  a  readily  adapted  scheme 
for  making  such  predictions,  and  to  give  (In  Appendix  A)  the  listing  of 
a  computer  program  for  calculating  the  stratification  parameter 
log1Q(us/h). 

The  advantage  of  using  a  hydrodynamic  model  to  compute  the  tidal  streams 
Is  that,  although  a  rough  estimate  of  the  parameter  u3/h  can  be  made  by 
using  tidal  charts,  this  method  cannot  be  used  in  regions  where  the 
tidal  currents  are  not  well  known.  However,  In  contrast  to  tidal  currents, 
tidal  elevations  have  been  extensively  studied  for  hundreds  of  years  and 
the  details,  for  almost  all  regions,  can  readily  be  found  in  the  open 
literature.  Taking  advantage  of  this  fact,  we  can  use  sea  level  Information 
as  Input  to  the  hydrodynamic  model,  l.e.  specify  the  rise  and  fall  of 
the  tide  around  the  boundary,  and  use  the  model  to  calculate  the  resulting 
Interior  tidal  currents  and  the  parameter  u*/h  In  regions  where  the 
tidal  currents  themselves  may  not  be  well  known. 


1  THE  HYDRODYNAMIC  EQUATIONS 

The  derivation  of  the  depth- integrated  form  of  the  two-dimensional 
hydrodynamic  equations  can  be  found  In  numerous  published  papers  and 
standard  texts  {[ a ,  9  s  10] ,  for  example).  However,  for  completeness,  a 
brief  (non-rigorous)  account  of  the  derivation  Is  given  here.  The  basic 
equations  are  (with  the  usual  notation): 
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where  right-handed  axes  have  been  chosen,  and  the  z-axls  Is  taken  as 
positive  downwards.  If  we  partition  each  velocity  component  Into  a  mean 
and  fluctuating  part  so  that 

u^  ■  uj  +  uj  , 

and  then  take  the  time  average  of  the  equations  over  an  Interval  that  Is 
long  compared  with  the  time  scale  of  the  fluctuations,  we  obtain 

if  *  5*7 (IIir)  *  5*7  (“1“r)  +  5*7 <5OT)  +  f=7  ■  - 1  Hr 

and 

^■  +  5^(5717)  +  377  (7177)  +  377  (7177)  - 


(N.B.  He  have  neglected  the  mean  quantity  non-linear  terms).  If  we  now 
drop  the  overbar,  and  use  eddy  coefficients  to  relate  the  turbulent 
stresses  to  the  gradients  of  the  mean  quantities,  then  we  obtain: 
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The  hydrostatic  equation.  Integrated  from  the  surface  (x3  *  -n)  down  to 
a  depth  x3  becomes 
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(we  have  neglected  the  horizontal  atmospheric  pressure  gradients). 
Therefore, 


Hr  =  9[Ht (X3+n)  p  +  (X3+n)  fer3 
=  gCHr^+  (Xs+n)  Hr 3* 

Hence, 


which  states  that  the  horizontal  pressure  gradient  has  two  parts:  one 
part  due  to  the  slope  of  the  free  surface,  and  the  other  part  due  to  the 
horizontal  variations  in  density. 

Hence,  the  pressure  term  on  the  right  hand  side  of  Eq.  1  becomes 
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[Eq.  3] 


The  next  stage  is  to  integrate  the  equations  vertically  and  express  the 
variables  in  terms  of  depth-mean  quantities.  For  simplicity  it  will  be 
assumed  that  the  water  is  well-mixed  vertically.  Under  these  circumstances 
p  =  p  (at  all  depths)  and  therefore 
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Consequently,  the  vertical  integration  of  the  pressure  terms  of  Eq.  3 
gives 


-a(d+n)  5n_  .  £ldM_  S£_ 

9^a  3Xi  p  2  3xj  * 

with  similar  terms  in  the  x2  momentum  equation.  To  integrate  the  velocity 
terms  over  the  depth  requires  the  use  of  kinematic  boundary  conditions 
at  the  top  and  bottom  boundaries. 
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If  4>(xi (X2 »x  3 )  ■  0  is  the  equation  of  the  surface  and  u^  =  (u1.U2.Uj) 
is  the  velocity,  then  the  surface  condition  is 

u  *  V*  =  -  §J  .  [Eq.  4] 

Since  the  surface  is  given  by 
x3  0  -n  (xi,x2)  , 

then 

<J>(xl,x2,x3)  =  n(xi,x2)  -  x3 
and 


Therefore,  Eq.  4  gives 

ft  +  Ul  1x7  +  U2  1x7  +  Us  =  0  at  X3  =  "n  •  [Eq*  5] 

At  the  bottom, 

x  3  =  d(xi,x2) 
and 

<J>(xi,x2,x3)  *  d(xi,x2)  -  x3. 

Therefore  ikV<|>  *  0  (i.e.  no  flow  through  the  bottom)  gives  that 

Ul  1x7  +  Uz  1x7  '  u3  =  0  at  x3  =  d.  [Eq.  6] 

If  the  velocity  terms  in  Eqs.  1  &  2  are  now  integrated  over  the  depth, 
and  the  kinematic  conditions  of  Eqs.  5  &  6  are  used,  then  the  momentum 
equations  become  (in  terms  of  depth-mean  quantities): 
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(Fy  and  Fv  are  the  components  of  the  surface  and 
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The  final  equation  that  is  required  is  the  depth-integrated  form  of  the 
continuity  equation,  which  will  be  used  to  compute  the  surface  elevation. 
If  v  is  an  arbitrary  volume  of  fluid,  fixed  in  space,  of  density  p  and 
surface  ares  A,  then 

L  I  pdv  =  -  j  pjj£- ds  +  /  Sdv  , 

1  v  A  v 

where  6  is  the  source  of  fluid/unit  volume/unit  time. 

This  gives 

/  [ft  +  dMpu)  -6]  dv  =  0, 
v 

and  hence 

+  div(pjj)  =  <5 . 


If  the  flow  is  assumed  to  be  incompressible,  this  reduces  to 
d1v(  u)  =  |  , 

i.e. 


3Ui  3u2  3u3 
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[Eq.  9] 


If  Eq.  9  is  now  integrated  vertically  over  the  water  column  and  use  is 
made  of  the  kinematic  conditions  of  Eqs.  5  &  6,  we  obtain 

It  +  7x7  Kd+ri>ui]  +  977  C(d+n)u23  -  (d+n)  |  .  [Eq.  io] 

This  equation,  along  with  the  momentum  equations  7  and  8,  are  the  ones 
solved  by  the  numerical  scheme.  In  the  present  analysis  the  horizontal 
density  gradient  terms  in  Eqs.  7  &  8  were  neglected. 
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FINITE  DIFFERENCES 


Equations  7,  8,  &  10  were  solved  using  centred  differences  on  a  regular 
grid  (Fig.  1).  (Note  that  the  grid  uses  left-handed  axes;  the  only 
effect  of  this  is  to  change  the  sign  of  the  coriolis  terms  in  Eqs.  7  &  8). 
Linear  interpolation  was  used  to  derive  the  values  of  variables  required 
at  points  other  than  grid  points.  A  leap-frog  scheme  was  used  for  the 
time  integration,  which  involves  the  storage  of  three  time  levels,  the 
time  derivatives  were  evaluated  using  levels  one  and  three,  while  the 
right-hand  side  of  Eqs.  7,  8  &  10  were  evaluated  at  the  middle  time 
level.  The  only  exceptions  were  the  diffusive  (friction)  terms;  these 
were  lagged  back  one  time  step  for  reasons  of  stability. 

As  an  example,  the  finite  difference  form  of  the  continuity  equation 
will  be  derived  (the  finite  difference  form  of  the  momentum  equations 
are  given  in  Appendix  A).  Let 
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be  written  as 
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then 

Ci  =  (d+n) 

where  6  is  the  source/unit  volume/unit  time. 

If  the  quantity  6'  is  the  total  source/grid  box/unit  time  then 

6'  =  (d+n)  (2Ax)  (2Ay)  6  . 

Therefore 
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p(d+n)  4AxAy  4pAxAy 

Consequently,  the  finite  difference  form  of  this  is  simply 

Ci  =  S'  (i  ,j)  /  [(4AxAy)  p(i,j)],  TL  =  2 

where  p(i,j)  denotes  the  value  of  p  at  the  grid  coordinate  (i,j) 
and  TL  =  2  signifies  time  level  2. 
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GRID  SCHEME  USED  IN 

THE  MODEL 
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Similarly, 
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This  is  the  equation  used  to  update  the  surface  elevations  (see  subroutine 
ETA3  in  Appendix  A);  comparable  expressions  can  be  derived  in  order  to 
update  the  U  and  V  velocity  components  (see  UVEL  and  VVEL  in  Appendix 
A).  In  certain  circumstances  (e.g.  when  a  point  lies  on  or  near  an  open 
boundary)  the  full  form  of  the  finite  difference  equations  are  not  used 
because  certain  terms  are  neglected.  The  points  were  coded  as  follows: 

For  elevation, 

(1)  outside  the  computational  grid,  do  not  update 
(i.e.  the  point  lies  Inland). 

(2)  a  boundary  value,  specify  through  a  boundary  condition. 

(3)  computational  point. 


For  velocity. 


(1) 

(2) 

(3) 

(4) 

(5) 


outside  the  grid  or  on  a  solid  boundary,  do  not  update, 
near  an  open  north-south  boundary,  nelgect  the  term 


involving  and  coriolis. 


near  an  open  east-west  boundary,  neglect  the  term 

£ 

in  rj-  and  coriolis. 

a  a  i 

both  of  the  above,  i.e.  near  an  open  corner  of  the 
grid. 

normal  interior  point,  update  using  the  full  form 
of  the  equations. 


The  finite  difference  form  of  the  equations  were  tested  by  applying  the 
model  to  simple  problems  having  known  solutions.  For  example,  the 
seiche  motion  of  a  rectangular  lake,  the  advance  of  a  wave  along  a 
rectangular  canal  (with  and  without  rotation),  the  wind-induced  set-up 
of  a  rectangular  bay,  and  similar  problems.  In  all  cases,  the  difference 
between  the  known  and  computed  solutions  did  not  exceed  a  few  percent  of 
the  analytical  solution.  The  application  of  the  model  to  the  interpreta¬ 
tion  of  some  field  observations  has  been  given  in  a  previous  memorandum 
(ill,  showing  a  good  agreement  between  the  observed  and  predicted  response 

In  the  present  study  the  amplitude  and  phase  of  the  surface  tide  is 
prescribed  around  the  open  boundary  of  the  model.  This  is  a  sufficient 
boundary  condition  to  drive  the  tidal  oscillations  (I.e.  tidal  currents 
and  elevations)  in  the  interior  of  the  model.  After  allowing  at  least 
five  tidal  cycles  for  the  system  to  approach  a  quasi -steady  state,  the 
subroutine  PARAM  (see  Appendix  A)  computes  the  value  of  log^uVh)  at 

each  elevation  point  by  averaging  the  surrounding  velocities.  At  each 
time  step,  and  for  each  grid  point,  the  present  value  of  the  parameter 
is  compared  with  the  previous  maximum  and  the  value  is  updated  if  a  new 
maximum  has  been  reached.  In  this  way,  the  maximum  value  of  log-.0(u3/h) 
is  calculated  at  every  computational  elevation  point.  u 


3  APPLICATION  TO  THE  SOUTHWESTERN  APPROACHES 
TO  THE  ENGLISH  CHANNEL 


As  an  example  of  the  scheme,  this  chapter  outlines  the  steps  Involved  in 
applying  the  model  and  calculating  the  stratification  parameter  in  the 
Southwestern  Approaches.  The  area  covered  by  the  model  is  shown  by  the 
insert  in  Fig.  2;  it  extends  from  the  Straits  of  Dover  in  the  east  and 
the  north  Channel  of  the  Irish  Sea  In  the  north  to  the  large  open  boundary 
south  of  Ireland  in  the  southwest.  The  eastern  and  northern  boundaries 
were  placed  at  narrow  sections  to  allow  a  better  definition  of  the 
boundary  conditions.  Open  boundaries,  like  the  one  to  the  southwest, 
should  be  avoided  whenever  possible  because  of  the  uncertainty  Involved 
in  specifying  the  tidal  constants  in  deep  water.  In  the  present  circum¬ 
stances,  however,  this  choice  of  the  boundary  was  unavoidable. 
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The  grid  points  of  the  model  and  the  solid  boundaries  (l.e.  the  coast- 
line)  are  shown  In  Fig.  3;  a  very  crude  grid  has  been  used  for  the 
purposes  of  demonstration,  to  obtain  more  accurate  results  It  would  be 
necessary  to  use  a  much  smaller  g~1d  spacing  (see,  for  example,  [7]). 

The  corresponding  water  depths,  defined  at  the  elevation  points,  are 
shown  In  Fig.  4.  The  amplitudes  and  phases  of  the  Mg  tide  around  the 

open  boundary  were  taken  from  Flather  [12J  and  are  shown  in  Fig.  5  for 
the  whole  of  the  model  region.  (The  boundary  values  are  given  in  the 
subroutine  BC).  The  major  uncertainty  (apart  from  the  effect  of  the 
crude  coastline  representation)  was  the  extrapolation  of  the  tidal 
boundary  conditions  along  the  large  open  boundary.  This  was  done  by  a 
process  of  trial  and  error  until  a  reasonable  agreement  between  the 
observed  and  computed  tidal  behaviour  was  obtained  for  the  western 
entrance  to  the  Channel  <  The  computed  phase  and  ampl 1 tude  distrfbuti ons 
are  shown  in  Fig.  6.  S'  fairly  good  agreement  was  achieved  in  the 
Approaches  to  the  Channel  and  in  the  Channel  itself,  and  also  within  the 
Bristol  Channel.  Note,  however,  that  the  agreement  between  the  observed 
and  computed  distributions  was  poor  throughout  most  of  the  Irish  Sea. 

The  amplitudes  shown  in  Figs.  5  and  6  are  for  the  Mg  component  alone.  As 

part  of  the  present  study  it  was  decided  to  investigate  the  variation  of 
the  stratif ication  parameter  (and  hence  the  predicted  frontal  location) 
during  the  spring/neap  cycle  of  the  tide.  The  proper  way  to  do  this 
would  be  to  specify  the  amplitude  of  both  Mg  and  Sg  around  the  open 

boundary,  to  run  the  model  to  simulate  about  one  month  of  elapsed  time, 
and  to  monitor  the  motion  of  the  frontal  location  during  the  spring/neap 
cycle.  A  simpler  way  is  to  run  the  model  for  three  cases: 

a.  Mean  tidal  conditions:  amplitudes  around  the  boundary 
given  by  Mg  alone. 

b  Spring  tides:  amplitudes  around  the  boundary  given 
by  (Mg  +  Sg). 

c.  Neap  tides:  amplitudes  around  the  boundary  given 
by  (Mg  -  Sg). 

Typical  ratios  of  the  MgiSg  amplitudes  for  the  region  were  obtained  from 

Heaps  [13],  and  the  open  boundary  amplitudes  were  scaled,  to  reproduce 
the  three  cases  given  above,  in  the  ratios  1.0  :  1.4  :  0.6. 

The  results  are  shown  in  Figs.  7,  8  and  9,  respectively.  They  suggest 
that  the  locations  in  which  the  front  might  possibly  form  cover  a 
considerable  proportion  of  the  area  at  the  western  entrance  to  the 
Channel.  However,  as  pointed  out  by  James  [6],  the  persistancy  of  a 
front  during  the  spring/neap  cycle  is  not  well  established.  It  is 
encouraging,  however,  that  the  calculated  mean  location  of  the  front 
(Fig.  7)  is  in  general  agreement  with  the  location  determined  from 
observations  [5,7,13]. 

Some  additional  tests  were  made  to  determine  the  likely  effect  of  storm¬ 
generated  currents  and  their  contribution  to  the  stratification  parameter. 
The  strongest  winds  in  this  region  are  directed  towards  the  east  [14] \ 
winds  of  Beaufort  force  8  or  greater  occurring  on  average  35-40  days/year. 
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FIG.  5 

THE  OBSERVED  AMPLITUDE  AND  PHASE  OF 
THE  M2  TIDAL  COMPONENT 

(taken  from  [12]) 


FIG.  6 

THE  COMPUTED  AMPLITUDE  AND  PHASE  OF 
THE  M2  TIDE 
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FIG.  11  CALCULATED  STRATIFICATION  PARAMETER  FOR  WEAK  TIDAL  CONDITIONS 
(boundary  amplitudes  equal  to  0.3  mean  tidal  amplitudes) 


FIG.  12  CALCULATED  STRATIFICATION  PARAMETER  FOR  WEAK  TIDES  PLUS 
THE  WIND  USED  IN  FIG.  10. 


This  corresponds  to  an  eastward  wind  of  around  16  m/s,  equivalent  to  a 
surface  stress  of  4  dyn/cm2.  To  simulate  the  effect  of  such  storms  the 
model  was  run  with  a  constant  eastward  wind  stress  and  the  wind-driven 
currents  computed.  As  In  the  previous  calculations  the  stratification 
parameter  1ogi0(u3/h)  was  estimated,  based  this  time  on  the  storm  driven 

flow.  (Note  that  this  calculation  Ignores  the  downward  vertical  mixing 
due  to  the  surface  stirring  of  the  wind.;  the  model  only  takes  account  of 
the  bottom-generated  turbulence).  The  wind-driven  flow  was,  almost 
everywhere,  too  weak  to  raise  the  stratification  parameter  to  the  critical 
value  (Fig.  10).  The  only  two  exceptions  were  a  location  off  the  south 
coast  of  Ireland  (where  there  was  an  easterly  flow  at  about  1  kn  (51  cm/s) 
in  the  shallow  water)  and  in  the  Channel  at  the  southern  entrance  to  the 
Irish  Sea  (where  the  flow  was  southward  at  about  1  kn).  Throughout  the 
whole  of  the  English  Channel  the  wind-driven  currents  alone  were  too 
weak  to  cause  the  generation  of  a  front.  In  general,  the  effect  of  wind 
driving  was  insignificant  in  comparison  with  the  tidal  effects,  even 
weak  tides  being  sufficient  to  mask  the  influence  of  the  wind.  To 
illustrate  this  point.  Fig.  11  shows  the  effect  of  weak  tides  (boundary 
amplitudes  set  to  0.3  of  the  mean  amplitudes),  while  Fig.  12  shows  the 
effect  of  weak  tides  plus  wind  on  the  stratification  parameter.  It  is 
clear  from  the  figures  that  the  tides  dominate  over  the  wind. 


CONCLUSIONS 


The  purpose  of  this  memorandum  has  been  to  present  a  readily  adaptable 
scheme  for  predicting  the  likely  locations  at  which  thermal  fronts  might 
be  formed  by  the  bottom  turbulence  associated  with  the  tides.  The  advantage 
of  the  method  is  that  it  requires  only  a  knowledge  of  the  behaviour  of 
the  tidal  elevations  around  the  open  boundary,  the  internal  dynamics 
within  the  region  of  interest  are  then  computed  numerically.  By  way  of 
an  example,  the  method  was  applied  to  an  area  that  includes  the  South¬ 
western  Approaches  to  the  English  Channel.  To  establish  the  computational 
grid,  digitize  the  depth,  and  code  the  grid  points  and  coastline  required 
about  five  working  days,  the  calculations  described  in  the  body  of  this 
memorandum  required  about  12  minutes  of  computer  time  each  on  a  Univac 
1106.  Consequently,  to  apply  the  scheme  to  a  new  region  and  make  all  of 
the  necessary  computations  would  require  about  10  man-days  of  effort 
(this  figure  should  be  maybe  doubled  or  trebled  if  the  individual  is  not 
familiar  with  the  model). 

In  the  present  study  the  model  was  adjusted  only  for  the  area  near  the 
western  entrance  to  the  English  Channel.  Better  results  would  probably 
be  obtained  if  a  higher  resolution  grid  was  used,  and  if  more  effort  was 
spent  on  calibrating  the  predicted  tide.  In  its  application  to  the 
Southwestern  Approaches  the  greatest  uncertainty  probably  arises  during 
the  extrapolation  of  the  tidal  amplitude  and  phase  along  the  open  boundary 
to  the  southwest. 
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VS  AR  (I  fJ  )=U. 

SB  AS  (I  fJ  »=0. 

CBARlI  fJ>  =  0. 

EB  Afi  (I  tJl=0. 

K1  (I  *J)=u. 

K2  (I  fJ  >=0. 

continue; 

NN=D 

NSUN=n 

RE  AO  IN  THE  DEPTHS  AT  S-POlNTS  <  H.  >• 

00  2  06  I=1»NX 

RE  AO  ( 5  »1  04  1  (  U<  It  JJ  »  J  =  l  »NT» 

FORMAT  (10F5.  0> 

READ  IN  THE  FRES.l  WATER  SOURCES  (H**3/S> 
00  208  1=1  .NX 

HE  AO  ( 5  f  1  Uft  >  CSIG(  If  J)fJ=l*NY) 


READ  IN  THE  WINU  SPEED  AND  DIRECTION  1H/S  AND  0E6REES* 
RE  AO  ( 5  fl  DO  1  UNO  /  CURRENT  CONVENTION) 

RE  AO  <  5  fl  DU  )  01  RN  ' 

COMPUTE  THE  COMPONENTS  OF  UlNO  STRESS. 

CC  =1.3*10.  **<-3.  » 

RQUA=1  J.  *10.  **  C-  3.  > 

01  RN  =  (  01  RN/180  .»  *3  .14159 
WX  =-  CC  *K  OWA*  IU  .*  *4  •*  WND*COSC  01  RN  )*  UND 
yy  =CC*  R0UA*1(%  **4.  *WMQ*SIN  (OIRAi)  *UNO 
01  RN  =<  01  RN  *1 80  .  >  /3  .1 4 1 59 


RE  AO  IN  THE  HO  HI  70  NT  AL  £00  Y  STRESSES. 

RE  AO  (5  ft  06  1  N 1  *N  2 
FORMAT  12F10.01 

RE  AO  <5  ft  01  1  IPLOTt 


171 

C 

I F  IPL0T1  •£ 0*  1  T<  OUTPUT  WILL  BE  PLOTTEO. 

172 

C 

173 

RE  AO  (5  *1 05  )  START 

17# 

105 

FORMAT !F  5*0) 

175 

C 

176 

C 

177 

C 

ESTABLISH  THE  INITIAL  CONDITIONS. 

178 

C 

179 

CALL  I  NIT 

180 

C 

181 

C 

OUTPUT  THE  INITIAL  VALUES. 

182 

C 

183 

C 

184 

UR  ITE?  6*  31) 

185 

31 

FORMAT (1  HI ) 

186 

32 

FORMAT  !/) 

187 

UR  ITE!6*  110) 

188 

1 10 

FORMAT  (5  X.14HINI  TIM.  VALUES///) 

189 

UR  IT E(  6«  111 )  OX 

190 

1 11 

FORMAT  (1UX*10H0E LT AX  r  *E10.3//) 

191 

UR  IT E(  6*  112)  OY 

192 

112 

FORMAT  (1  OX  .10HUELT  AY  =  tElO.J//) 

193 

UR  IT  E!  6*  113)  OT 

194 

113 

FORMATUOXvlOHOELTAT  =  »El0.3////> 

195 

UR  ITE!  6*  114)  TOR 

196 

1 14 

FORMAT  (1  OX* 25H 80  TT  CM  FRICTION*  TOR  =  *E10.3//> 

197 

if)  IT  E!  6*  115)  RRLAT 

198 

115 

FORMAT  <1  OX  »12HLATI  1UDE  =  *F5.0//) 

199 

UR  IT E!  6*  11 6)  UNO  *0  JR N 

2  JO 

116 

FORMAT  (1  OX  *14HUI  NO  SPEED  =  *F  5.  0*  10X*  18HUIN0  DIRECTION 

2  01 

UR  ITE!  6*  11?)  N1»N2  /  •  F5.0/) 

FORMAT  (1  OX  *32H  HORIZONTAL  EDDY  STRESSES  Nl'r) 

2  02 

117 

203 

4E10.3//)  /  tE10.3*10X*6)4i  2  =  *1 

14?  ITE C 6*  118)  RK1  *RK2  • 

204 

205 

118 

FORMAT  <1  OX  .37HH0RI  ZONTAL  EDDY  01  FFUSIVITIES  K1  =  • 

206 

*E1  0.  3.  10X.6HK2  =  £10.3//) 

207 

WRITE?  6*31) 

208 

C 

289 

C 

OUTPUT  THE  UATER  DEPTHS. 

210 

C 

211 

UR  ITE (6*  120) 

21 2 

120 

FORMAT  11  OX  *21HUA  TER  DEPTH  IN  METRES*//) 

213 

C 

214 

DO  215  I=i*NX 

215 

WRITE*  6*32) 

216 

215 

UR  ITE!  4*  119)  <  0<  I*  J)  *J=1*NY> 

217 

119 

FORMAT  !1X*16F7  .1  ) 

218 

WRITE! 6*  31) 

219 

C 

220 

UR  ITE?  6*  121) 

221 

121 

FORMAT  !10X»28HFRESH  UATER  SOURCES  IN**3/S)  *//) 

222 

C 

223 

00  216  1=1. NX 

224 

UR  ITE!  6*  32) 

225 

216 

UR  ITE!  6*  119)  ISIG!  I*J)*J=1*NY) 

226 

C 

227 

C 

24 


228 

C 

229 

C 

CONVERT  TO  CSS  UNITS. 

230 

c 

231 

00  207  I=1**X 

232 

00  207  J=l#NV 

233 

0<  I.  J)  =0  (1  *J>*10.«  *2  • 

234 

H<  I*  J»  1)=0  (I  tJ  > 

235 

H<  I*  J.  2>=D  (I  *J» 

236 

SIG<  £»  J)  =SIG  (I  *J  >*10.* *6. 

237 

207 

CONTINUE 

238 

C 

239 

C 

240 

C 

241 

C 

24  2 

C 

243 

C 

244 

CALL  OUTPUT 

245 

C 

246 

c 

2  47 

c 

*****  START  THE  MAIN  COMPUTATIONAL 

LOOPS  **** 

248 

c 

**  **  * 

**** 

249 

c 

**  •*  * 

*••* 

250 

c 

***** 

•  **• 

251 

c 

**  **  * 

**** 

252 

c 

**  **  * 

**** 

253 

c 

254 

TO=TOR 

2  55 

THAT  =0  0.  3 

256 

c 

257 

00  3  00  NN=1*MST£P 

258 

c 

259 

t=floa  T<  NN  >  *  OT 

260 

TI  M=  T/  36  UU  . 

261 

TC=T  If*  /l  2.4  224 

262 

c 

263 

c 

264 

TO  R=  TO  *T  MA  X/EX  P(  TC  > 

265 

c 

266 

c 

267 

c 

268 

CALL  8  C 

269 

c 

270 

CALL  ETA 3 

271 

c 

272 

c 

273 

CALL  UVEL 

274 

c 

275 

CALL  VVEL 

276 

c 

2  77 

CALL  PAR  AM 

278 

c 

279 

c 

280 

c 

*****  all  variables  have  now  been 

UPOATEO**** 

281 

c 

232 

c 

UPDATE  THE  TIME  INUEX. 

233 

c 

284 

00  720  I=1*NX 

25 


285 

00  72C  J=ltUY 

286 

U(  If  Jf  1»=U([  *J  *2  ) 

287 

U(  I»  Jt  2)  =UCI  *J  f 3  ) 

288 

VC  If  Jt  i)=VCI  fj  f2  ) 

289 

VC  If  Jf  2)=V(I  tj  fJ  ) 

290 

ET  AC  If  Jf  1)=ETA  «  tj  #2  ) 

291 

ETAC  If  Jf  2)=£TA  Cl  f  J  f3  » 

292 

HC  If  Jt  1)=H  Cl  f  J  f 2  > 

293 

H<  If  Jf  2)  =H  Cl  f  J  f3  » 

294 

720 

CONTINUE 

295 

C 

296 

C 

297 

c 

TEST  FOR  OUTPUT. 

298 

c 

299 

IFCTC  JLE.  START!  GO  TO 

300 

300 

c 

3  01 

IF  (MOD  (N.Nf  NANS!  .ME.  0! 

GO  TO  699 

302 

c 

303 

c 

PR  IN  T  OCJ  T  RE  SU  LT  S. 

304 

c 

305 

CAUL  OUTPUT 

306 

699 

CONTINUE 

307 

C 

308 

C 

309 

3  CO 

CONTINUE 

310 

c 

311 

IF  Cl  PL  OT  1  .EQ.  V  t  GO  TO 

1255 

312 

CALL  PLO  TCO.  fO.f  9991 

313 

125  5 

CO  iNTINUE 

314 

C 

315 

C 

316 

STOP 

317 

END 

1 

SU  BR  OU  n  ME  I  NI  T 

2 

C 

3 

C 

SET  THE  INITIAL  CONDITIONS. 

4 

C 

5 

C 

6 

C 

INSERT  THE  INITIAL 

VALUES. 

7 

C 

8 

00  100  I  =  lfNX 

9 

DO  100  J=lfNY 

10 

DO  100  K=l»3 

11 

IF  Cl  SC  If  J)  .EQ  .  1> 

GO  TO  100 

12 

ROWCIf  Jf  K)  =  l  . 

13 

1  00 

CONTINUE 

14 

c 

15 

RETURN 

16 

END 

26 


SUBROUTINE  OUTPUT 


50 


C 


c 

c 

OUTPUT  THE  RESULTS. 

c 

URITEI6.  ill 

31 

FORNAT  tl  HI  1 

32 

FORMAT  1/1 

C 

C 

OUTPUT  THE  TIME  IN  HOURS. 

t 

UR  IT  E4  6.  25)  TI  R»  TC 

25 

FORMAT  Cl  OX  *25H£L  APSE  0  TINE  IN  HOURS 
%  26  HN  UN  BE  R  OF  TIOAL  CYCLES  =  .F*.2/) 

=  .F6.2.50X* 

C 

C 

UR  IT  E  4  6,100) 

100 

FORMAT  I5X.18HSURFA IE  ELEVATIONS/) 

00  1  50  I=1.NX 

UR  IT  E  4  6*  32) 

150 

UR  IT  E4  6.  101 )  4  ETA(  I*  J.2)  *J=1  *NY) 

101 

/* 

FORNAT  C1X*16F7.1  ) 

L 

c 

URITC46.il) 

URITE4  6*  102) 

102 

FORMAT  HOX.31HO»COIf>ONENT  HORIZONTAL 

VELOCITY.///) 

00  151  I=1«NX1 

WRITE (6,  32) 

151 

r 

URITE4  6.  lul)  4  U4  I*  J*  2)  t  J=1  «NY) 

V 

UR  IT  E4  6.  31 ) 

UR  IT E 4  6.  103) 

103 

r 

FORNAT  UQX.31HV-COf*»ONENT  HORIZONTAL 

VELOCITY./// 1 

00  1  52  I  =1  .NX 

UR  IT  E  4  6.  32  ) 

152 

C 

C 

UR  IT  £4  6.  101)  4V4  I*  >2).J=1.NY1> 

c+*+ 

URITE4  6.31) 

URITE4  6*  104) 

10« 

FO  RN  AT  a  ox  .  24H  ST  RA  Q  FI  CA  TI  ON  PAR  AH  ETER  »///  ) 

V 

00  1  53  1=1. NX 

URITE4  6.  32) 

153 

UR  IT  E4  6.  101)  I  SI  I.  J»  2)  »J=1  .NY) 

C 

C 

♦♦♦♦♦♦♦♦♦♦♦♦  ♦.  ♦♦♦♦♦♦  +»*»♦♦♦ 

C 

t 

IF  II PL  OT 1  .EQ.  0)  60  TO  125 

CALL  GRAPH 

125 

CONTINUE 

27 


I 


H 


L  '' 


T 


57 

C 

58 

C 

59 

c 

+■+■++■+■+■*++++  +  ++++  ♦+♦+  ♦♦♦♦♦♦  ♦♦  ♦♦  ♦♦  ♦ 

• 

60 

c 

61 

RE  TURN 

62 

ENO 

- 

t 

SUBROUTINE  tiC 

2 

C 

J 

C 

SPECIF  r  THE  TIDAL  AHPLITUQES  A  NO 

4 

C 

PHASES  ArtOUNO  THE  OPEN  B0UN0ART. 

5 

C 

6 

C 

r 

SC  ALE-  Q.  60 

8 

a 

9 

ETA(  1.  6*  31=2  00.*  SC  AL  E*  SIN  1 0.  5058  *T  El -0  . 

10 

<r 

11 

c 

12 

00  V=  SI  NC  0.  5058  *T  IN  1 

13 

c 

14 

c 

15 

ET  A(  9*  12  *31=250.  *SCALE*00V 

16 

ET  AC  10  *12*31=300  •*  SC  ALE*  00  V 

17 

ET  AC  tl  *1  2*  31  =350  •*  SCALE*D0  V 

ia 

c 

19 

c 

20 

A=  120. 

21 

TT=0.5Q5B*TIH 

22 

c 

23 

ET  AC  8,  1*  3)  =A*SCALE  *SlNITT-0.  0175*150*1 

24 

ET  AC  9,  1,  3)  =A  *S  CA  LE  *S  IN  C  T  T-0.  01 75 *140.1 

25 

ET  AC  10  *1  *3  1=A*SCALE*SINCTT-0.0175*130.) 

26 

ET  A<  11  ,1  ,3  >  =  A*  SC  ALE*  S  INC  TT -0-0 17  5*  125.1 

27 

ET  AC  12  *1  *3  1=A* SCALE* SI NCTT -0.0175*  120.) 

28 

ETAC  13  *1  *3>  =  A*  SC  ALE+SINCTT -0.0175*  120.) 

29 

ETAC  14  *1  *3  >  =  A*  SC  ALE*SINCTT  -0.0175*115.) 

30 

ETAC  15  *1  *3  1= A* SC AL£*SINCTT -0.0 175*110.) 

31 

ET  AC  16  *1  *3  )=  A*  SC  AL  E*  SI  N(  TT  -0  .0 17  5*  10 5.  ) 

32 

ETAC  17  *1  *3  )=A*  SC  ALE*  SI  NCTT  -0.0175*  100.) 

33 

ETAC  17  *2  *3  )=  A*  SC  ALE*  SI  NCTT  -0.0 17  5*  100.) 

34 

ETAC  17  *3  *3  >= A*  SC  ALE*  SI  NCTT -0.0 17 5*  100.) 

35 

ET  AC  17  *4  *3  )=  A*  SC  ALE*  SI  NCTT  -0.0 17  5*  100.) 

36 

ETAC  17  *5  *3)=A«SCALE*SINCTT-0.017S*100.) 

37 

ETAC  17  *6  *3  )=  A*  SC  AL  E*  SI  N<  TT  -0  -0175*  100. 1 

38 

ETAC  17  *7  *3  )= A*  SC  AL  E*  SI  NC  TT  -0  .0 17  5*  100.  > 

39 

ET  AC  17  *8  *3  )=  A*  SC  AL  E*  SINC  TT  -0  .0175*  100.  ) 

40 

c 

41 

RE  TURN 

42 

ENO 

28 


SUBROUTINE  t  TA3 


* 

I 


K  . 


UPOATE  THE  SURFACE  ELEVATIONS. 


00  1  00  I  =  i»NX 
00  100  J=1*NV 

IF  (I  E(  I*  J>  .EQ.  1)  GO  TO  100 
IF  (I  £t  It  J)  .EQ.  2)  GO  TO  101 

COMPUTATIONAL  POINT. 

C1=SIG(I  tJ>/(4.«DX*0Y*R0UCIt  Jt  2)  > 

C2=-  CC  H(  It  Jt2)  *U  Cl  *1  •  J*2  >)  *CJ  (I  »J  *2  >■*  CH  tl  ~1  *J  »2I*HC  I«  J*2)  ) 

*  *U  Cl  -1  *J  *2)  )  /( A.  *0X> 

C3=-  <<  H(  It  J*2J*H  Cl  »J*1  t2l)  *V  <1  »J  t2H  CH  Cl  tJ-1  .2)  *H<  It  Jt2>> 

*  *V  (I  tJ  .2  I)  /<4.  «0Y) 


ETAC  I»  J.  J)=ETA  (I  ,J  tit*  2.  *0  T*  (C  1*C2*C3) 
CONTINUE 

H(  It  Jt  3)  =ETA  Cl  tJ  *3  HOC  It  J) 

CO  NTINUE 

RETURN 

END 


SUBROUTINE  UVEL 


1 

2  C 

J  c  UPDATE  The  u  components. 

4  c 

5  REAL  Mat  #*62  *Mrt3  *Md4  »MB5*Hd6  .MB7 

6  00  ICO  I:t»NXl 

7  00  IPO  J=1*NY 

a  c 

9  IF  Cl  Ut  I.  J>  .EQ  •  1)  i>0  TO  100 

10  C 

11  C  OUTSIDE  THE  GR  IU  OR  ON  *  BUY. 

12  C 

13  we  j=  ux 

14  C 

15  *  6=  -G  *<  H(  I.  J.2)  *H  II  ♦It  0*2  i)  *<  ETAt  1*1*  J.2)  ~€  T A  ( I  *  J*2)  )/ <4.  *0 X) 

16  MB7=-G*<  (H(I  »J  •  ?)♦*  1*1 «J* 2)  >**2.>  *(  ROW  (  I*  1 .  J  *2  > -ROW  « I  *  J  *2  >  I 

17  S  *(  ROW(  It  o»2)  ♦ROW  (I  *i  »J»2)  )  **  t-  1.  >/ <8.*0X> 

18  C 

19  C  NEAR  AH  OPEN  £  -W  JCY  OR  A  CORNER. 

20  C 

21  IF(ILMI'J)  .£0.  5  .OR.  IU(I.J)  .  EQ  .  4  )  GO  TO  20 

22  C 

23  C  NORMAL  EQUATION 

24  MB  1=*N  1*  IH  (I  ♦!  *0  *1  >*  <U<I*l  *J  *1  >-  U(  I*  J*1  )  l-H<  I*  J*l)*IU«I*  J*  1) 

25  $  -U  (I  -l  *0*11)  1/  (  4  •*  OX  **  2.  1 

26  GO  TO  21 

27  C 

2b  2C  Mb  1=0. 

29  21  CONTINUE 

30  C 

31  C 

32  C  NEAR  A, \l  OPEN  N  -S  8  UY  OR  A  CORNER. 

33  C 

34  IF(IUfl.J)  .EQ.  2  .OR.  IU<I»J>  .  EQ  .  4 )  GO  TO  22 

35  C 

36  C  NORMAL  POINT. 

37  «)  shi=o. 25* <  h<  i.  j.  i)  ♦hu  »j*i  *t  >♦«<  i*i*  j*i*  i)»HCin»Jtin 

38  H)  SH  ?=  0.  23*f  HI  I*  J-  1*11  ♦«  <1  »J  *1  )♦  HI  I*l»  J*  1)  ♦«  f  I  *1  *  J— 1  *1)  ) 

39  MB  2=  *N  2*  (HUSHl  *(  U<  I*  J+l»  1>  -UII  *0.1  I)  -HDSH2*<  U<  1*  J*  1 1 -t)  «X  *  J-l  *11) 

40  S  )/  (4  .*  OY  **2.  » 

41  C 

42  GO  TO  23 

43  C 

44  22  MB  2=0. 

45  23  CONTINUE 

46  C 

47  C 

48  C 

49  C  NEAR  AN  OPEN  N-S  8UY  OR  A  CORNER. 

50  C 

51  IF  (I  U<  I*  J)  .£0.  2  *0  R«  IUfI*J>  .EQ.  4)  60  TO  24 

52  C 

53  U1  =U(I  .0  *1 1 

54  U2=0.2  5*  <V  (I  *0-1  *1  >♦  V<  I*  J*  1)  *V  (I  *1  *J *1  »*V(  1*1. 0-1*1)  > 

55  «B4:-0  .5  •(  ROW(  I*  J*  1)  ♦ROW  (1*1  *J  *1  ))  *TOR*U  1*SQRT  I Ul**2  .♦U3**2.  ) 

56  GO  TO  25 


30 


K0)0)0t<4*4»is*4'<l<<l'4'>4a|0'^0  6-6-^60-0-frUIUi 

UNHOOCt'JO'UHUMH'  ^CB'JffUlfloiNHOOa 


m 


57  C 

24  M8  4=-0  «S  *(  R0U(  I#  J.  II  ♦R0W(I*1  *J*1  >1  •TOR* 

SU<  If  J*  1)  *ABS  (U  (I  fj  til) 

25  CONTINUE 
C 
C 

C  NEAR  AN  OPEN  N-S  SOY  OR  A  CORNER. 

C 

IF  (I  U(  I*  J)  ,EQ.  2  .0  R.  IU(IfJ)  .  EQ.  4>  60  TO  26 
C 

MB5=*F  *(  H(  It  Jf  2)  *M  (I  *lf  Jf2>)  •(  V<  It  Jt 2)  ♦Vtl-M  .  Jt2> 

«*V  <1  *1  fj-l  f 2  )*  V(  If  J-  It  21  >/«. 

C 

SO  TO  27 

26  MB  5=0. 

C 

2  7  CONTINUE 

C 

c 

U<  If  jf  31  =  (0.5*<H  (I  tJ  fl)*H<  I*lf  jf  1)  »*U<If  J.l) 

*♦2.*  OT  *(  MB  1  ♦  MB  2*  NO  3*  M B  4*  MB  5*  MB  6*  MB  7)  1/ *0  .S*<  HC I  f  Jf  3) 
t+H  (1*1  .J  f3>>  > 

C 
C 

100  CONTINUE 

RE  TU  RN 
END 


31 


SUBROUTINE  WE  L 


1 

2 

C 

3 

C 

A 

C 

5 

6 

7 

8 

C 

9 

C 

10 

c 

11 

12 

c 

13 

14 

c 

15 

16 

c 

17 

18 

19 

c 

20 

c 

21 

c 

22 

c 

23 

24 

25 

26 

27 

28 

c 

29 

30 

c 

31 

20 

32 

c 

33 

21 

34 

c 

35 

c 

36 

C 

37 

38 

c 

39 

40 

41 

c 

42 

43 

c 

44 

22 

45 

c 

46 

23 

47 

C 

48 

c 

49 

c 

50 

51 

c 

52 

53 

54 

c 

55 

56 

UPDATE  THE  V  COMPONENTS. 

REAL  Kdl  tM82  tHB3  tH  B4  tMB5tHB6  tHB7 
DO  100  I=1»NX 
DO  ICG  J=ltNYl 

OUTSIDE  THE  GR 10  OR  ON  A  BOY. 

IF  (I  Vdt  J)  .  EQ.  1)  GO  TO  100 
HB3=WY 

MS  6=  -G  ♦<  H(  I.  J*l»  2)  «H  (I  *J»2  >)  *<  ETA<  I*  J*  1*2)  -C  TA  f  I  *  J  t2>  I /<4.  *0  Y» 

MB  7=  -G  *«  H<  I.  Jt  2>  *H  (I  t  J-*l  t2  >>  *« 2.  •<  ROUdt  J*lt2)-R0U(I  »J*2)  ) 

* /(  ROUdt  Jt2»  ♦ROU  (I  tJ*lt2>>  /(  8.*DY) 

AE  AR  AN  OPEN  E-U  BUY  OR  A  CORNER. 

IF  (I  V<  It  J)  .EQ.  3  .OR.  IVdtJI  .EQ.  A)  60  TO  20 
K)SHl=u.25*(H(  I.  Jt  11  ^HII  tjtl  tl  >♦«<  I*  It  JMt  1)+Htl*-1  tJtll » 

M3  SH2=ll.  25*<H(  It  It  Ji  1J+MCI  -1  tJ*l  tl  >*HC  It  J*  it  1  »♦«(  I  tJ  1 1>  ) 

KB  1-  1*  (H0SH1  *<  7(  1+  It  Jt  1)  -VCI  tJ  tl  ))  ~H0SH2*t  Vdtjt  D-VCI-l  tJ  till  » 

*  ✓<  4.  *0  X*  *2.1 

60  TO  21 

H8  1=0. 

CO  NT  I  Nl£ 

NEAR  AN  OPEN  N-S  HOY  OR  A  CORNER. 

IF  (IV(I*J)  .EQ.  2  .OR.  IVdtJI  .EQ.  41  GO  TO  22 

KB2=*N2*  (H  (I  tJ*l  tl  )*(V(ItJtltU-V(  It  Jt  1  I  >-H<  I  tJ  tl )  • 

$  (V  (I  tJ  tl  >-V<It  J- It  1)  >)/<«.  *0Y*  *2. ) 

GO  TO  23 

MB  2=0. 

CONTINUE 

NEAR  AN  OPEN  E-U  BOY  OR  A  CORNER. 

IF  tl  V  f  1 1 J I  .EQ.  3  .OR.  IVdtJI  .EQ*  4)  60  TO  20 

U1  =0  .2  5*  CU  1 1  -1  tJtl  BU(I-lt  Jt-ltll  ♦U(I.J*ltl)*UUtJtlH 
U2=V(I  ,Jtl) 

f%*4=  -0  .5  *(ROU(  It  Jt  ll  ♦ROUCI  tjtl  tl  >1  *T0R*U2»SQRT tUl**2 .♦U2*«2.  » 

60  TO  25 


57 

C 

58 

24 

MB  4=  *0  .5 

59 

25 

CONTINUE 

60 

C 

61 

c 

62 

c 

NEAR  AN 

63 

c 

64 

IF  (I  VI  I* 

65 

C 

66 

MB  5=  *f  *4 

67 

t«U  (I  »J«1 

68 

C 

69 

60  TO  27 

70 

26 

*5=0. 

71 

C 

72 

27 

continue 

73 

C 

14 

C 

75 

C 

76 

VI  If  Of  31 

77 

S  (MB1  *M  B2 

78 

c 

79 

c 

80 

100 

CONTINUE 

81 

c 

82 

RETURN 

83 

END 

NEAR  AN  OPEN  E-U  BOT  OR  A  CORNER  • 

HFlIVlIvJk  .EQ.  3  JDR.  IVlIfUl  .EQ.  41  60  TO  20 


1 

2  C 

3  C 
A  C 

5  C 

6 

7 

8  C 

9 

10  C 

11 

12 

13  C 

14 

15  C 

16 
17 

10  C 

19 

20  C 

21  100 

22  C 

23 

24 


SUBROUTINE  PAR  AN 

COMPUTE  THE  STRATIFICATION  PARAMETER* 

00  1  00  I  -2  tNXl 
00  1  00  J=2  *NY1 

IF  II  SI  I*  J)  .EQ.  II  60  TO  100 

C<  If  J*  2>=I0.5*«JII  -1  »d.2>*U(I.  Jf2l  >**2. 

t*0J*<vUf>ii2)«va>J»2n«*2.) 

Cl  I.  J»  2)  =SQRTI  C«  I.  J.  2)  » 

ST  RA  T=CC  I*  J*2)  «3.  /D«I  fJ  1*0.  001 
Cl  I*  Jt  2)  =AL0G1  Of  ST  RATI 

F(C(I  fJf2l.6T.SII  «J  *2  II  SllfJ  f2l=CIIf  J*2I 

CONTINUE 

RETURN 

END 


SU  BR  OU  TI  NE  6  Ra  Ph 


1 

2  C 

3  c  oo  the  plotting  on  a  cal comp *96u. 

4  c 

5  C 

6  C 


7  C 

OLEV 

THE  CONTOUR  SPACING 

8  C 

HGTC 

ThE  HEIGHT  OF  THE  CONTOUR  LABEL  IN  INCHES. 

9  C 

INOC 

PERMITS  CONTOUR  LABELS  TO  VARY  DURING  A  SINGLE  PLOT. 

10  C 

LA  dC  .GE  .  0 

NO.  OF  01  GITS  IN  CONTOUR  LABEL  DECIMAL  PART. 

11  C 

=  -  i 

NO  DECIMAL  PART. 

12  C 

=  -3 

NO  LABEL  PRINTED  ON  THE  CONTOUR  LINE. 

13  C 

LW6T=1 

CONTOUR  DRAWN  8Y  ORDINARY  LINE. 

19  C 

-2 

HEAVY  LINE. 

15  C 

=3 

□  ASHE  0  LINE. 

16  C 

NA  ac=i 

VARIES  1-10  CAN  BE  USED  TO  SMOOTH  THE  CONTOURS. 

17  C 

NLEV 

SUPPLIED  BY  GETLEV.  I  T  IS  THE  NO.  OF  CONTOURS  ORAUN. 

18  C 

NX  *N  Y 

U0RKIN3  DIMENSIONS  OF  THE  ARRAY  BEING  CONTOUREO. 

19  C 

Xi  PL  .v  IP  L 

COORDINATES  OF  THE  7  ORIGIN.  I.  E  Z(l.l). 

20  C 

XL  PL  »Y  LP  L 

COORDINATES  OF  Z<NX*NY>. 

21  C 

Z 

ARRAY  TO  BE  CONTOURED. 

22  C 

ZLEV 

ARRAY  SUPPLIED  BY  GETLEV. 

23  C 

29  C 

25  C 

THE  F0LL0UING  ARE  THE  PLOT  OPTIONS. 

26  C 

27  C 

IP  LO  Tt  =0 

N  0  PL  OT  . 

28  C 

IP  LOTI  =1 

PLOT  ELEVATIONS  AND  VELOCITY  VECTORS  ONLY. 

29  C 

IP  LOTI  =2 

ALSO  PLOT  SALINITY. 

30  C 

31  C 

32  C 

33 

00  X=0.  79  07 

39 

00  Y=  0.  99  5  9 

35  C 

36 

FCX=2.  *0  OX  *F  LO  AT  (N  XI  ) 

37 

FC  Yz2.  *0  OY*F  LO  AT  <N  Y1  > 

38 

X^HFT=2.  *1)  0  X  *F  LO  AT  (NXl+2. 

39 

YSHFT=2.  *OUY  .FLOAT  <NY  J-.2. 

90  C 

91  C 

SET  THE  PLOTTING  PARAMETERS. 

92  C 

93 

INOCzl 

94 

LWGTzl 

45 

HGTC=0  .4 

46 

LA  BC  =-  1 

47 

NARC  =  l 

48  C 

49  C 

50  C 

IF  NN=0  THEN 

CONTOUR  THE  BOTTOM  TOPOGRAPHY. 

51  C 

52 

IF  <NN  JL  0.  0)  GO  TO  100 

53  C 

59 

GO  TO  101 

55  C 

56  100 

CONTINUE 

■  - 

m>  4K>  * 

34 


57  C 

58  C 

59  C 

60 
61 

62  C 

63 

64 

65 

66 

67  1 02 

68  C 

69  C 

70 

71  C 

72  C 

73  C 

74 

75  C 

76  C 

77  C 

78 

79 
SO 

81  C 

82  C 

83  C 
34 
85  C 
36 

87 

88 

89  C 

90 

91 

92  C 

93 

94 

95 

96 

97 

98 

99 

100  C 

101 

102  C 

103  101 

104  C 
1  15 
106 

107 

108 

109  C 

110 

111  C 

112  C 

113  C  _ 


CONTOUR  THE  BOTTOM  DEPTHS. 

CALL  N£UPLT<4.  »S  AA  S  3»0> 

CALL  F  ACTOR<  0.  5) 

00  102  I  =  1*NX 
00  102  J=1*NY 
Z(  I*  J)  =0  (I  »  J  )  /1C  .*  *2  . 

IFIIS(I.J)  .£0.  1)  2(ItJ)  =  10.**35. 
CONTINUE 

SET  CONTOUR  SPACING  AT  5U0  M  • 
0L£V=50C. 


CALL  G  £T  l£ V <  Z *  NX  »N  Y*  DLEV  »ZLEV»NLEV  > 

CONTOUR  THE  UE  PT  H. 

CA  LL  PL0T<4.  ,7.. -3  » 

CALL  CONTUR(Z*NX  »M  Y»  0.  tO.t  FCX»FCY*  ZL EV »L ABC* LUGT *NL£V t HGTC * 
*NARC  ) 

OR  AO  THE  BOUNDARY. 

CALL  BOUNO 

CALL  SYM80L(d.  »- 2.  «u.28*16HDX  OY  DT.0..*16> 

CALL  S  YMHQL ( 5. 0*  — 2  •#  0.  28 *23HT0 K  UX  WY  F» 

4  0..  »*2  3  ) 

POX=OX*l  0.  **<-5.  > 

PO  Y=CY  *1  U.  *•  (-5.  > 

CA  lL  N  UM  BE R <  U.  »-3.  .0  .28  »  PD  X*  0.  *«-2» 

CALL  NUM8ER(2.*-3.  «0.28»PDY»0.»»2) 

CALL  N  UM  8ER<4.  *-3.  »0.28t0T*0.*  *2 1 
CALL  NUNBER<6.  »- 3.  «O.28tT0R»  0.r*4l 
CALL  NUM8ER<8.  •- 3.  tO  .28*UX  »0  .t  *2  > 

CALL  NUMBERUQ  .•  -3  .»0.  28»UY»0.  r*2J 
CA  LL  N UM BE R <  1 2  .*  -J  .♦  0.  28*F  *0  .»+6> 

CALL  PLOTtO.  »0.»9991 

CONTINUE 

IF<IPL0T1  .EQ.  1>  CALL  NE  WPL  T<  4t  »SWA  •  *  3  *0» 

IF(IPL0T1  .EQ.  2  .OR.  IPLOTi  .EQ.  3) 

♦  CALL  N  EW  PLT  (  3#  *SUA*»3*Q> 

CA  LL  FACT0R(O.5> 

CALL  PL0T<4.  »7.t  -3  1 

CONTOUR  THE  SURFACE  ELEVATIONS. 
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114  00  110  I=1.NX 

115  00  110  J=1.NY 

116  21  I.  J»=ETA<I  *J*2  > 

117  IF(IE<  I«J)  .E0.  1)  Z  ( I  •  J  >  =  10  .*  *3  5. 

118  110  CONTINUE 

119  C 

120  C  SET  THE  CONTOUR  SPACING  AT  100  CM. 

121  C 

122  OL  EV=1 00  • 

123  C 

124  CALL  6ETLEVCZ.NX  *N  Y*OLEV»ZL£V*NL£V> 

125  C 

126  CALL  CONTUR(Z*NX  ,N  Y.  0.  .0  . »  FC  X.  FC  Y*  ZLE V  * L  ABC*  LUGT  .NLEV.H6TC  * 

127  *NARC> 

128  C 

129  C 

130  C  DRAW  TtC  VELOCITY  VECTORS. 

131  C 

132  USCAL£=0  .02 

133  C 

134  00  1  20  1=1*  NX1 

135  DO  120  J=1  *NY1 

136  C 

137  UN  00  =U  (I  *  J  *2  )*J*2  •  +  IN  I  *  «!♦  1*2)  **  2.  ♦  V<  I*  J.2)  *  *  2.  ♦  V<  I  ♦  1*  J,  2)  **2. 

138  C 

139  IFtUHOO  .EQ.  0.)  GO  TO  120 

140  C 
14  1  C 

142  XC  0=  DO  X*2.  *F LO  AT  (I  -1  >*DDX 

143  >C0=00Y*2.  *FLOAT  (J -1  )*00Y 

144  C 

145  CALL  SYHBOLCXCO*  YC  0.  0.  16.4*0  .*“11 

146  C 

147  UB  AR=0  .5*<Ut  I*  J.  21  *0(1  *J*1  .2  1) 

148  VBAR=0.5*<  VII.J.2H  Vfltltj»2)i 

149  C 

150  C 

151  XD  SH=XCO*UBAR*USCAl£ 

152  YDSH=YCO-*VBAR*USCAl£ 

153  C 

154  CALL  PLOT<  XOSH  *YOSH*  2) 

155  C 

156  120  CONTINUE 

157  C 

158  C 

159  C  PLOT  T>€  BOUNDARY. 

160  C 

161  CALL  BOUND 

162  C 

1  o3  CALL  SVHB0L<1.  r- 2.  «0.28»4HTC  *0  .»  -Ml  1 

164  CALL  NUNBERll.  r-  3.  <0  .28*  TC.O  .* -»2> 

165  CALL  S  YN0OL<3.  r- 2.  «0.28*5HTIN  *0.  r*51 

166  CALL  NUHBER13.  *-3.  «0.28*Tin*0. 

167  C 

168  C  CONTOUR  THE  STRATIFICATION  PARAMETER. 

169  C 

170  IF  (<  IP  LOTI  .EO.  21  .OR.  (IPL0T1  .EQ.  3>1  60  TO  130 


171  60  TO  131 

172  1  30  CONTINUE 

173  C 

17*  CA  LL  PLOTCO.  #YSHFT  r*3> 

175  C 

176  00  132  1=1  »NX 

177  00  132  J=1'NY 

178  ZC  I*  J)  =S  Cl  *J  '2  ) 

179  132  CONTINUE 

180  C 

Idle  SET  CONTOUR  SP  AC  IN  6  AT  1 

182  C 

183  0LEV=1  JO 
18*  C 

185  CALL  GET  LEV  C  Z*  NX  »N  Y*  DLEV  'ZLEVtNLEV  > 

186  C 

187  C 

188  CALL  C  ON  TURC  Z*  NX  *N  Y*  0*  *0  •  »  FC  X*  FC  Y»  2LEV  tL ABC*  LUGT  *NLEV*H6TC  t 

189  SNARC) 

190  C 

191  C 

192  131  CONTINUE 

193  C 

19*  C  DRAW  T»€  BOUNDARY. 

195  C 

196  IF  1C  IP  LO  Tl  .EQ.  2)  .OR.  (IPL0T1  *E  0.  3>>  GO  TO  150 

197  GO  TO  151 

198  150  CONTINUE 

199  C 

200  CALL  8  01)  NO 

201  C 

202  151  CONTINUE 

203  C 

20*  CALL  P  LO  T  (  0.  tO  •  1 99  91 

205|C 

206|  C 

207  RETURN 

208  END 
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SUBROUTINE  BOUND 

OKAUS  THE  COASTAL  BOUNDARY. 


CO M  ON  /C  iO /  XA  (2 91  *Y  Al  20  >•  XB  (20)  *YB(  20) •  XC (201  * YCC  20>* 

*»(2  0>  »Y0(20>.  XE  (2  0)  *Y£(  20  It  XF  (20)  »YF(20> 

DA  TA  ( X  A(  J )  *J=l  *7  )/U.  .2  .A  .2  .A  .8  .7  »B  .7  *0 »•  1. / 

DA  U  (V  A(  Jl  *J=1  *7  )/ 3.  *3.«  1.  *1  .*0.  *0  .*  1./ 

DA  TA  ( X B(  J)  *J=1  *9  i/Qm  .2. 4. 2.4  *8.7  •«  .7*10.3*  10.3.0..  I./ 

DA  TA(YB(  Jl  *J=1  *9)79.  *9. *7.  *7  .*5.  •  5  ••  0.  *0. 1 1.  / 

OA  TA  (X  C(  J)  »J=1  *2  0)  A).  *0.8*  0.8*  2.  4*  2.  4*  5.5*  5. 5*7.1  • 

*  7.  1*  8.7*  8. 7*  10.3*1  Ok  3*  11.9*1 1.9*  13.4»13.4*15.0*0.*1.7 
OA  TA  ( Y  C<  J)  *J  =  l  *2  3)  /1 1.  *11.  *13.  .13.  *15.  *15.  *11.  *11.  *13.  *13.  » 
*9.  *9.*  13. *1 3.*  11  .*  11  .*9.  *  9  •*  0.  *1./ 

DATA(XD<J>  *J=1  *15)  A>.*0.B*0.8»2.4*2.4*7.1*7.1*11.9*11.9* 

*13  .4  *1  3.  9*15.0  *15.  0*  0.  *1./ 

OA  TA  (YOf  J)  *J=1  *1  3)  /I  7.  *17.  *19.  *19.  *21.  *21.  • 

* 23  .*  23  .»  21  •*  21  •»  13  •*  13  •* 9.  »0.*1./ 

OATA  (XE(  J)  *J=1  *15)  A3.  4*  13  .4  *1 6.  6*  16  .6  *18.  2*  18.  2* 

*16  .6  *1  6.  6*  19.8  *19.  8*  22.9*22.9*  25.3  *0.*1./ 

OA  TA  (Y  E<  J)  ,J=1  *15)  Ai  3.  9*  23  •*  23  •«  21  •*  21  •*  1 9  •* 

*19.*  17  .*  17.*  11  •*  11  .*  1S.*15.*0.  *1  •/ 


DRAW  THE  BOUNOARY  OF  THE  GRID. 


CALL  PLOTCO.  *0.*3) 

CALL  PL0TC25.3  *0.*2) 
CALL  PLO  7(25.30*23.9*2) 
CALL  PLO  TIO.  *23.  90  *2  ) 
CALL  PLOTIO.  *0.*  2) 


CALL  L  INEiXA  *YA*5*  1*  0*0) 
CALL  L  INE(XB*YB*7»  1*  0*0) 
CALL  L  INE(XC*YC«  18  *1  *0*0) 
CALL  LINE(X0.Y0*  13*1  *0*0) 
CALL  LINE(XE.YE*  13  *1*0*0) 
RE  TURN 


uiuiuiuiyiwiui#*^«#*****uwwwwwwwwwNrvjNNwwrjMNrvHHHHHHPHHH 

»ui»WMH34a>JO'Ui*uruHO'Oag^O'Ui*UMHO-ca<j»ui»uNH04ai^»ui«UNM04a^»'Vii«uNitf 


27.78 
33.34 
74.5344 
0017 
0013 
00  2  2  00 
0300 

ii  ii  13  ii  mu 
n  u  1551  n  m 

11 11  55  5511111 
11 11  55  11 11 11 1 
11115511  11111 
111151  11  11111 
1115  51  11  11111 
2555  5511  11111 
25555111  11121 
25  55  51 15  55  521 
2555  55  5551511 
2555  55  5551111 
255555  11 11111 
25555511 11111 
25555511 11111 
43  3333  3311111 
1111 11  11  1111 
11 11  15  11  11 11 
1111  55  51  1111 
111155  51  1111 
11 11  51  11  11 11 
1111  55  11  1111 
11151111  1111 
25555511  1111 
2555  51  11  1111 
2555  11 15  5551 
2555  5555  5551 
25  55  55  551111 
2555  55  551111 
2555  51 11  1111 
255551  11  1111 
25555551  1111 
43  33  33  31  till 
11 11 12  Ulllll 
1111133111111 
11113  3  3311111 
111133  33  11111 
11 113311 11111 
111133  3111111 
111331  Ulllll 
23333331 11111 
2333  3311  11121 
23333113  33  321 
2333  33  33  33  321 
2  333  33  3331311 
23333333  31111 
233333  11 11 11 1 
23  3  3  3311 11111 
2  3  3  3  33  3311111 
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57  222222  2211111 

58  11 11 12  11  11 11 1 

59  1111131111111 

60  111133  3311111 

61  1111333311111 

62  11 1133  U  11  111 

63  111133  3111111 

69  11 1331  11  11  111 

65  23  33  33  3111111 

66  2333  33  11  11 121 

67  2333  31  13  33  321 

68  23333333  33  321 

69  2333  33  33  31 31 1 

70  2  3  33  33  33  31111 

71  23  3  3  3311  11111 

72  2333331111111 

73  2  3  33  33  3311111 

74  22  22  22  2211111 


75 

• 

• 

• 

• 

18. 

• 

• 

• 

• 

76 

• 

• 

n 

• 

• 

• 

• 

64. 

36. 

• 

• 

78 

• 

• 

79 

• 

• 

• 

9  . 

54. 

36. 

7. 

• 

60 

• 

• 

61 

• 

• 

• 

18  • 

53. 

27. 

5. 

• 

82 

• 

• 

83 

• 

• 

• 

7  . 

73. 

• 

• 

• 

84 

• 

• 

85 

• 

• 

• 

36. 

73. 

7. 

• 

• 

86 

• 

• 

87 

• 

• 

36. 

9. 

• 

• 

• 

• 

88 

• 

• 

89 

36. 

18. 

73. 

82. 

100. 

45. 

16. 

• 

• 

90 

• 

• 

91 

137. 

11  7. 

95  . 

95. 

92. 

60. 

• 

• 

• 

92 

• 

9. 

• 

93 

146. 

91. 

73. 

110. 

82. 

• 

• 

36. 

18. 

27. 

94 

45. 

1  1. 

• 

95 

120. 

10  4. 

1 10. 

110. 

44. 

55. 

73. 

66. 

82. 

62. 

96 

36. 

7. 

• 

97 

128. 

14  6. 

120  . 

126. 

1  10  . 

92. 

84. 

55. 

18. 

98 

11. 

• 

99 

146. 

16  0. 

120  . 

14  6. 

1  20  . 

106. 

88. 

18. 

9. 

100 

• 

• 

101 

920. 

16  0. 

150. 

110. 

1  20. 

92. 

• 

• 

• 

102 

• 

• 

103 

24  70.183  0. 

3  60  . 

126. 

1  06  . 

92. 

• 

• 

• 

104 

• 

• 

• 

105 

4350.330  0. 

18  50  . 

360. 

1  37. 

130. 

92. 

11. 

• 

106 

• 

• 

• 

107 
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